Reads were first de-multiplexed with ea-utils (fastq-multx) with an allowed barcode mismatch of 1 (Default). De-multiplexed reads were then subjected to an alignment against the UniVec vector contamination database using the adapter removal tool within ea-utils (fastq-mcf). Barcodes were extracted by remove the last 12 nucleotides from each read.
Next reads were subjected to an alignment against the Oligo library using a DNA aligner. We used Bowtie2…
We are going to first import the merged counts tables that is outputted from the Nextflow pipeline. This table already has the raw counts summarised from the alignments and merged into a single matrix.
# Import counts
cd8.counts = readRDS('../data/CD8-counts.rda') %>%
rownames_to_column("Id") %>%
mutate(Id = recode(Id, oligo_1899 = "MART1_ELA_a",
oligo_4281 = "MART1_ELA_b",
oligo_913 = "MART1_a",
oligo_3295 = "MART1_b",
oligo_1887 = "CDK4mut_a",
oligo_4269 = "CDK4mut_b",
oligo_1888 = "CDK4wt_a",
oligo_4270 = "CDK4wt_b",
oligo_1889 = "RPS3Amut_a",
oligo_4271 = "RPS3Amut_b",
oligo_1890 = "RPS3Awt_a",
oligo_4272 = "RPS3Awt_b",
oligo_1895 = "MANSC1mut_a",
oligo_4277 = "MANSC1mut_b",
oligo_1896 = "MANSC1wt_a",
oligo_4278 = "MANSC1wt_b")) %>%
column_to_rownames("Id")
# Convert to long format
cd8.counts.long = cd8.counts %>%
rownames_to_column("Id") %>%
pivot_longer(cols = -c('Id'), names_to = 'sample', values_to = 'counts')
# Import CD4 metadata
cd8.metadata = readRDS('../data/CD8-metadata.rda')
# Import MultiQC results
cd8.report = readRDS('../data/CD8-multiqc.rda')
# Import oligo annotations
oligo.genes = readRDS('../data/oligo-genes.rda')
# Mean
mean(cd8.report$raw_total_sequences)
## [1] 14162045
# Median
median(cd8.report$raw_total_sequences)
## [1] 13417612
# Range
range(cd8.report$raw_total_sequences)
## [1] 8434046 23626015
cd8.report %>%
select(Sample, reads_unmapped, reads_mapped) %>%
mutate(Sample = str_replace(Sample, '-bamstats', '')) %>%
pivot_longer(cols = c("reads_unmapped", "reads_mapped"), names_to = "feature", values_to = "reads") %>%
ggplot(aes(x = Sample, y = reads, fill = feature)) +
geom_col() +
coord_flip() +
theme_light(base_size = 11) +
theme(legend.position = "top") +
scale_fill_aaas() +
xlab("") +
ylab("Number of reads")
cd8.report %>%
select(Sample, reads_unmapped_percent, reads_mapped_percent) %>%
mutate(Sample = str_replace(Sample, '-bamstats', '')) %>%
pivot_longer(cols = c("reads_unmapped_percent", "reads_mapped_percent"), names_to = "feature", values_to = "reads") %>%
mutate(reads = reads / 100) %>%
ggplot(aes(x = Sample, y = reads, fill = feature)) +
geom_col() +
coord_flip() +
theme_light(base_size = 11) +
theme(legend.position = "top") +
scale_fill_aaas() +
scale_y_continuous(labels = scales::percent) +
xlab("") +
ylab("Percent of reads")
cd8.counts %>%
ggplot(., aes(x = log10(mock_a), y = log10(mock_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("Mock replicate #1 (log10)") +
ylab("Mock replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
cd8.counts %>%
ggplot(aes(x = log10(id3_a), y = log10(id3_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("1D3 replicate #1 (log10)") +
ylab("1D3 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
cd8.counts %>%
ggplot(aes(x = log10(cdk453_a), y = log10(cdk453_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("CDK4/53 replicate #1 (log10)") +
ylab("CDK4/53 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
cd8.counts %>%
ggplot(aes(x = log10(d1_10_a), y = log10(d1_10_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("1:10 replicate #1 (log10)") +
ylab("1:10 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
cd8.counts %>%
ggplot(aes(x = log10(d1_100_a), y = log10(d1_100_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("1:100 replicate #1 (log10)") +
ylab("1:100 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
cd8.counts %>%
ggplot(aes(x = log10(d1_300_a), y = log10(d1_300_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("1:300 replicate #1 (log10)") +
ylab("1:300 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
cd8.counts %>%
ggplot(aes(x = log10(d1_1000_a), y = log10(d1_1000_b))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("1:1000 replicate #1 (log10)") +
ylab("1:1000 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..))
# Calculate the variance between oligo-replicates
oligo.variance = cd8.counts %>%
rownames_to_column("Id") %>%
mutate_if(is.numeric, funs(log10)) %>%
left_join(oligo.genes, by = "Id") %>%
group_by(oligoGroup) %>%
summarise(mock_a = diff(mock_a),
mock_b = diff(mock_b),
cdk453_a = diff(cdk453_a),
cdk453_b = diff(cdk453_b),
id3_a = diff(id3_a),
id3_b = diff(id3_b),
d1_10_a = diff(d1_10_a),
d1_10_b = diff(d1_10_b),
d1_100_a = diff(d1_100_a),
d1_100_b = diff(d1_100_b),
d1_300_a = diff(d1_300_a),
d1_300_b = diff(d1_300_b),
d1_1000_a = diff(d1_1000_a),
d1_1000_b = diff(d1_1000_b))
## `summarise()` has grouped output by 'oligoGroup'. You can override using the `.groups` argument.
# Plot histogram of variance across samples
oligo.variance %>%
pivot_longer(-oligoGroup, names_to = "Sample", values_to = "difference") %>%
ggplot(aes(x = difference)) +
geom_histogram(bins = 100) +
theme_bw() +
facet_wrap(.~Sample) +
xlab("Difference (log10)")
# Match the sample ID's between samples
cd8.counts = cd8.counts[,match(cd8.metadata$Sample, colnames(cd8.counts))]
# Make DEseq2 object
ddsMat.cd8 <- DESeqDataSetFromMatrix(countData = cd8.counts,
colData = cd8.metadata,
design = ~Group)
## converting counts to integer mode
# Get normalize counts
ddsMat.cd8 <- estimateSizeFactors(ddsMat.cd8)
ddsMat.cd8 <- estimateDispersions(ddsMat.cd8, fitType='local')
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
# Get normalized counts
cd8.norm <- counts(ddsMat.cd8, normalized = TRUE) %>%
as.data.frame() %>%
rownames_to_column("Id")
# Write tables to disk
write_tsv(cd8.counts, "tables/cd8-counts.tsv")
write_tsv(cd8.norm, "tables/cd8-counts-norm.tsv")
# Create a long table
cd8.norm.long = cd8.norm %>%
gather(sample, norm.counts, -Id)
# Pre-normalized
p1 = cd8.counts.long %>%
ggplot(aes(x = sample, y = log10(counts + 1))) +
geom_boxplot() +
theme_light(base_size = 11) +
theme(axis.text.x = element_text(angle = 90)) +
xlab("") +
ylab("Counts (log10)") +
ggtitle('Un-normalized data')
# After normalized
p2 = cd8.norm.long %>%
ggplot(aes(x = sample, y = log10(norm.counts + 1))) +
geom_boxplot() +
theme_light(base_size = 11) +
theme(axis.text.x = element_text(angle = 90)) +
xlab("") +
ylab("Normalized counts (log10)") +
ggtitle('Median-normalized data')
# Plot side-by side
p1 + p2
cd8.norm.long %>%
ggplot(aes(x = log10(norm.counts))) +
geom_histogram(aes(y=..density..), binwidth=0.25, colour="black", fill="white") +
geom_density(alpha = 0.2, fill = "#FF6666") +
facet_wrap(sample~.) +
theme_light(base_size = 11) +
ylab("Density") +
xlab("Normalized reads per oligo (log10)")
# Get average between the replicates
id3.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
ID3_mean = rowMeans(select(., id3_a, id3_b))) %>%
mutate(M = log2(ID3_mean / Mock_mean),
A = (log2(ID3_mean) + log2(Mock_mean)) / 2)
# Plot scatterplot
p1 = id3.mean %>%
ggplot(aes(x = log10(Mock_mean), y = log10(ID3_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("Mock normalized counts (log10)") +
ylab("1D3 normalized counts (log10)")
# Plot MA plot
p2 = id3.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " 1D3/Mock)")))
# Plot side-by side
p1 + p2
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("mock_a", "mock_b", "id3_a", "id3_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "Mock")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
# Get results (ID3 / Mock)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4762 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 0, 0%
## LFC < 0 (down) : 5, 0.1%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
id3.mock.res = results(ddsMat.subset, alpha = 0.05)
# View table
id3.mock.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05, abs(log2FoldChange) > 1)
## Id baseMean log2FoldChange lfcSE stat pvalue
## 1 MART1_ELA_a 1616.221 -2.592509 0.1668882 -15.53441 2.029294e-54
## 2 MART1_b 2359.897 -1.175441 0.1121242 -10.48339 1.029893e-25
## 3 MART1_ELA_b 1750.376 -4.113238 0.1543848 -26.64277 2.170495e-156
## padj
## 1 4.831750e-51
## 2 1.634784e-22
## 3 1.033590e-152
# Show quick plot of significance ()
DESeq2::plotMA(id3.mock.res, ylim = c(-3, 3))
# Get average between the replicates
cdk4.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
CDK4_mean = rowMeans(select(., cdk453_a, cdk453_b))) %>%
mutate(M = log2(CDK4_mean / Mock_mean),
A = (log2(CDK4_mean) + log2(Mock_mean)) / 2)
# Plot scatterplot
p1 = cdk4.mean %>%
ggplot(aes(x = log10(Mock_mean), y = log10(CDK4_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("Mock normalized counts (log10)") +
ylab("CDK4 normalized counts (log10)")
# Plot MA plot
p2 = cdk4.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " CDK4/Mock)")))
p1 + p2
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("mock_a", "mock_b", "cdk453_a", "cdk453_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "Mock")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
# Get results (CDK4 / Mock)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4759 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 3, 0.063%
## LFC < 0 (down) : 4, 0.084%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
cdk4.mock.res = results(ddsMat.subset, alpha = 0.05)
# View table
cdk4.mock.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05 & abs(log2FoldChange) > 1)
## Id baseMean log2FoldChange lfcSE stat pvalue
## 1 CDK4mut_a 772.9214 -2.890157 0.2133074 -13.54925 8.003909e-42
## 2 CDK4mut_b 620.7388 -4.422590 0.3231989 -13.68380 1.268798e-42
## padj
## 1 1.904530e-38
## 2 6.038211e-39
# Show quick plot of significance ()
DESeq2::plotMA(cdk4.mock.res, ylim = c(-5,5))
Oligo Id’s for the positive control dropout’s are: - oligo_1899 (MART1_ELA) - oligo_4281 (MART1_ELA) - oligo_912 (MART1) - oligo_3294 (MART1)
# Get average between the replicates
id3.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
CDK4_mean = rowMeans(select(., cdk453_a, cdk453_b)),
ID3_mean = rowMeans(select(., id3_a, id3_b))) %>%
mutate(M = log2(ID3_mean / CDK4_mean),
A = (log2(ID3_mean) + log2(CDK4_mean)) / 2)
# Plot scatterplot
p1 = id3.mean %>%
ggplot(aes(x = log10(CDK4_mean), y = log10(ID3_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("CDK4 normalized counts (log10)") +
ylab("1D3 normalized counts (log10)")
# Plot MA plot
p2 = id3.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " 1D3/CDK4)")))
p1 + p2
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("cdk453_a", "cdk453_b", "id3_a", "id3_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "CDK4")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
# Get results (ID3 / CDK4)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4750 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 4, 0.084%
## LFC < 0 (down) : 5, 0.11%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
id3.cdk4.res = results(ddsMat.subset, alpha = 0.05)
# View table
id3.cdk4.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05 & abs(log2FoldChange) > 1)
## Id baseMean log2FoldChange lfcSE stat pvalue
## 1 CDK4mut_a 804.7891 2.955708 0.1968698 15.01352 5.988381e-51
## 2 MART1_ELA_a 1499.7890 -2.466632 0.1630949 -15.12390 1.126636e-51
## 3 MART1_b 2439.4898 -1.244110 0.1100652 -11.30339 1.262446e-29
## 4 CDK4mut_b 730.6316 4.667656 0.2748360 16.98342 1.089508e-64
## 5 MART1_ELA_b 1793.0922 -4.150197 0.1537456 -26.99393 1.741353e-160
## padj
## 1 7.111203e-48
## 2 1.783841e-48
## 3 1.199324e-26
## 4 2.587582e-61
## 5 8.271426e-157
# Show quick plot of significance ()
DESeq2::plotMA(id3.cdk4.res)
# Get average between the replicates
d10.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D10_mean = rowMeans(select(., d1_10_a, d1_10_b))) %>%
mutate(M = log2(D10_mean / Mock_mean),
A = (log2(D10_mean) + log2(Mock_mean)) / 2)
# Plot scatterplot
d10.scatter = d10.mean %>%
ggplot(aes(x = log10(Mock_mean), y = log10(D10_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("Mock normalized counts (log10)") +
ylab("D10 normalized counts (log10)")
# Plot MA plot
d10.ma = d10.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 5) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " D10/Mock)")))
d10.scatter + d10.ma
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("mock_a", "mock_b", "d1_10_a", "d1_10_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "Mock")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
# Get results (D10 / Mock)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4759 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 0, 0%
## LFC < 0 (down) : 7, 0.15%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
d10.mock.res = results(ddsMat.subset, alpha = 0.05)
# View table
d10.mock.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05 & abs(log2FoldChange) > 1)
## Id baseMean log2FoldChange lfcSE stat pvalue
## 1 CDK4mut_a 800.8005 -2.507273 0.2217680 -11.30584 1.227801e-29
## 2 MART1_ELA_a 1621.4895 -2.555543 0.1588604 -16.08672 3.161408e-58
## 3 CDK4mut_b 624.8460 -4.227006 0.2828515 -14.94426 1.697680e-50
## 4 MART1_ELA_b 1735.4737 -4.350605 0.1653703 -26.30826 1.542788e-152
## padj
## 1 1.460776e-26
## 2 7.522572e-55
## 3 2.693086e-47
## 4 7.342126e-149
# Show quick plot of significance ()
DESeq2::plotMA(d10.mock.res)
# Get average between the replicates
d100.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D100_mean = rowMeans(select(., d1_100_a, d1_100_b))) %>%
mutate(M = log2(D100_mean / Mock_mean),
A = (log2(D100_mean) + log2(Mock_mean)) / 2)
# Plot scatterplot
d100.scatter = d100.mean %>%
ggplot(aes(x = log10(Mock_mean), y = log10(D100_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("Mock normalized counts (log10)") +
ylab("D100 normalized counts (log10)")
# Plot MA plot
d100.ma = d100.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 5) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " D100/Mock)")))
d100.scatter + d100.ma
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("mock_a", "mock_b", "d1_100_a", "d1_100_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "Mock")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.
## final dispersion estimates
## fitting model and testing
# Get results (D100 / Mock)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4763 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 9, 0.19%
## LFC < 0 (down) : 13, 0.27%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
d100.mock.res = results(ddsMat.subset, alpha = 0.05)
# View table
d100.mock.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05 & abs(log2FoldChange) > 1)
## Id baseMean log2FoldChange lfcSE stat pvalue
## 1 CDK4mut_a 837.8370 -2.121543 0.1884076 -11.260387 2.058525e-29
## 2 MART1_ELA_a 1685.6399 -2.207934 0.1242659 -17.767816 1.254810e-70
## 3 CDK4mut_b 716.7704 -2.263236 0.5349687 -4.230595 2.330739e-05
## 4 MART1_ELA_b 1916.8408 -2.656286 0.1128993 -23.527917 2.113133e-122
## padj
## 1 3.268252e-26
## 2 2.988329e-67
## 3 1.264367e-02
## 4 1.006485e-118
# Show quick plot of significance ()
DESeq2::plotMA(d100.mock.res)
# Get average between the replicates
d300.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D300_mean = rowMeans(select(., d1_300_a, d1_300_b))) %>%
mutate(M = log2(D300_mean / Mock_mean),
A = (log2(D300_mean) + log2(Mock_mean)) / 2)
# Plot scatterplot
d300.scatter = d300.mean %>%
ggplot(aes(x = log10(Mock_mean), y = log10(D300_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("Mock normalized counts (log10)") +
ylab("D300 normalized counts (log10)")
# Plot MA plot
d300.ma = d300.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 5) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " D300/Mock)")))
d300.scatter + d300.ma
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("mock_a", "mock_b", "d1_300_a", "d1_300_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "Mock")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.
## final dispersion estimates
## fitting model and testing
# Get results (D300 / Mock)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4762 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 8, 0.17%
## LFC < 0 (down) : 7, 0.15%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
d300.mock.res = results(ddsMat.subset, alpha = 0.05)
# View table
d300.mock.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05 & abs(log2FoldChange) > 1)
## Id baseMean log2FoldChange lfcSE stat pvalue
## 1 MART1_ELA_a 2044.892 -1.071998 0.1144178 -9.369152 7.311744e-21
## 2 MART1_ELA_b 2355.140 -1.239435 0.1006218 -12.317753 7.269179e-35
## padj
## 1 1.740926e-17
## 2 3.461583e-31
# Show quick plot of significance ()
DESeq2::plotMA(d300.mock.res)
# Get average between the replicates
d1000.mean = cd8.norm %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D1000_mean = rowMeans(select(., d1_1000_a, d1_1000_b))) %>%
mutate(M = log2(D1000_mean / Mock_mean),
A = (log2(D1000_mean) + log2(Mock_mean)) / 2)
# Plot scatterplot
d1000.scatter = d1000.mean %>%
ggplot(aes(x = log10(Mock_mean), y = log10(D1000_mean))) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab("Mock normalized counts (log10)") +
ylab("D1000 normalized counts (log10)")
# Plot MA plot
d1000.ma = d1000.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 5) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " D1000/Mock)")))
d1000.scatter + d1000.ma
# Create new DESeq2 object with poor samples removed
ddsMat.subset = ddsMat.cd8[,which(colnames(assay(ddsMat.cd8)) %in% c("mock_a", "mock_b", "d1_1000_a", "d1_1000_b"))]
ddsMat.subset$Group = droplevels(ddsMat.subset$Group)
ddsMat.subset$Group = relevel(ddsMat.subset$Group, "Mock")
# Run stats
ddsMat.subset <- DESeq(ddsMat.subset)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
# Get results (D1000 / Mock)
results(ddsMat.subset, alpha = 0.05) %>%
summary()
##
## out of 4762 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 0, 0%
## LFC < 0 (down) : 2, 0.042%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Get table of results
d1000.mock.res = results(ddsMat.subset, alpha = 0.05)
# View table
d1000.mock.res %>%
as.data.frame() %>%
rownames_to_column("Id") %>%
filter(padj < 0.05 & abs(log2FoldChange) > 1)
## [1] Id baseMean log2FoldChange lfcSE stat
## [6] pvalue padj
## <0 rows> (or 0-length row.names)
# Show quick plot of significance ()
DESeq2::plotMA(d1000.mock.res)
# Scatter plot
d10.scatter + d100.scatter + d300.scatter + d1000.scatter
# MA plot
d10.ma + d100.ma + d300.ma + d1000.ma
# Rank
cd8.norm.long %>%
group_by(Id) %>%
summarise(sum = sum(norm.counts),
mean = mean(norm.counts),
median = median(norm.counts),
sd = sd(norm.counts),
geomean = gm_mean(norm.counts),
cv = sd(norm.counts) / mean(norm.counts)) %>%
mutate(rank = rank(mean)) %>%
ggplot(aes(x = rank, y = log10(mean))) +
geom_point() +
theme_minimal()
# Get mean counts per oligo
cd8.mean = cd8.norm.long %>%
group_by(Id) %>%
summarise(sum = sum(norm.counts),
mean = mean(norm.counts),
median = median(norm.counts),
sd = sd(norm.counts),
geomean = gm_mean(norm.counts),
cv = sd(norm.counts) / mean(norm.counts))
# Get quantiles
quantile(cd8.mean$median , seq(from = 0, to = 0.10, by = 0.005))
## 0% 0.5% 1% 1.5% 2% 2.5% 3%
## 0.000000 0.318799 1.062006 1.978373 9.955168 27.003683 74.638207
## 3.5% 4% 4.5% 5% 5.5% 6% 6.5%
## 142.062296 218.769373 319.850709 396.491177 445.067499 500.165253 551.465274
## 7% 7.5% 8% 8.5% 9% 9.5% 10%
## 592.921535 651.214172 689.947011 741.235201 767.365712 801.337192 831.117392
# Histogram with 5% cutoff
cd8.norm.long %>%
group_by(Id) %>%
summarise(total.count = sum(norm.counts),
mean.count = mean(norm.counts),
median = median(norm.counts),
geomean = gm_mean(norm.counts),
sd.count = sd(norm.counts)) %>%
ggplot(aes(x = log10(median))) +
geom_histogram(binwidth = 0.25, colour = "black", fill = "white") +
geom_density(alpha = 0.2, fill = "#FF6666") +
geom_vline(xintercept = log10(396.491177), color = "red") +
theme_light(base_size = 11) +
xlab("Normalized reads per oligo (log10)")
# Get ID's within the limits
## Cutoff = 5%
idx <- cd8.norm.long %>%
group_by(Id) %>%
summarise(sum = sum(norm.counts),
mean = mean(norm.counts),
median = median(norm.counts),
geomean = gm_mean(norm.counts),
sd = sd(norm.counts) ,
c.v = sd/mean) %>%
filter(median >= 396.491177) %>%
pull(Id)
# Filter table
cd8.norm.fil = cd8.norm %>%
filter(Id %in% idx)
# Write filter table to disk
write_tsv(cd8.norm, "tables/cd8-counts-norm-filtered.tsv")
# Get dimensions
dim(cd8.norm)
## [1] 4764 15
dim(cd8.norm.fil)
## [1] 4525 15
# Create long dataframe
cd8.norm.fil.long = cd8.norm.fil %>%
gather(sample, norm.counts, -Id)
# Plot boxplot
cd8.norm.fil.long %>%
ggplot(aes(x = sample, y = log10(norm.counts + 1))) +
geom_boxplot() +
theme_light(base_size = 11) +
theme(axis.text.x = element_text(angle = 90)) +
xlab("") +
ylab("Normalized counts (log10)") +
ggtitle('Median-normalized data w/ filtering')
# Histogram all
cd8.norm.fil.long %>%
ggplot(aes(x = log10(norm.counts))) +
geom_histogram(aes(y=..density..), binwidth=0.25, colour="black", fill="white") +
geom_density(alpha = 0.2, fill = "#FF6666") +
facet_wrap(sample~.) +
theme_light(base_size = 11) +
ylab("Density") +
xlab("(Filtered) Normalized reads per oligo (log10)")
# Get average between the replicates
id3.mean = cd8.norm.fil %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
CDK4_mean = rowMeans(select(., cdk453_a, cdk453_b)),
ID3_mean = rowMeans(select(., id3_a, id3_b))) %>%
mutate(M = log2(ID3_mean / CDK4_mean),
A = (log2(ID3_mean) + log2(CDK4_mean)) / 2)
# Plot MA plot
id3.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab(expression(paste("Fold change ", '(', log[2], " 1D3/CDK4)")))
### Titration data
# 1:10
d10.mean = cd8.norm.fil %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D10_mean = rowMeans(select(., d1_10_a, d1_10_b))) %>%
mutate(M = log2(D10_mean / Mock_mean),
A = (log2(D10_mean) + log2(Mock_mean)) / 2)
d10.plot = d10.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 1) +
xlim(7, 14) +
xlab("") +
ylab("")
d10.plot
# 1:100
d100.mean = cd8.norm.fil %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D100_mean = rowMeans(select(., d1_100_a, d1_100_b))) %>%
mutate(M = log2(D100_mean / Mock_mean),
A = (log2(D100_mean) + log2(Mock_mean)) / 2)
d100.plot = d100.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 1) +
xlim(7, 14) +
xlab("") +
ylab("")
d100.plot
# 1:300
d300.mean = cd8.norm.fil %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D300_mean = rowMeans(select(., d1_300_a, d1_300_b))) %>%
mutate(M = log2(D300_mean / Mock_mean),
A = (log2(D300_mean) + log2(Mock_mean)) / 2)
d300.plot = d300.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 1) +
xlim(7, 14) +
xlab("") +
ylab(expression(paste("Fold change ", '(', log[2], " D300/Mock)")))
d300.plot
# 1:1000
d1000.mean = cd8.norm.fil %>%
mutate(Mock_mean = rowMeans(select(., mock_a, mock_b)),
D1000_mean = rowMeans(select(., d1_1000_a, d1_1000_b))) %>%
mutate(M = log2(D1000_mean / Mock_mean),
A = (log2(D1000_mean) + log2(Mock_mean)) / 2)
d1000.plot = d1000.mean %>%
ggplot(aes(x = A, y = M)) +
geom_point(colour = alpha("grey", 0.5)) +
geom_point(colour = alpha("red", 0.7), data = . %>% filter(Id %in% c("MART1_ELA_a", "MART1_ELA_b"))) +
geom_point(colour = alpha("lightblue", 0.7), data = . %>% filter(Id %in% c("MART1_a", "MART1_b"))) +
geom_point(colour = alpha("orange", 0.7), data = . %>% filter(Id %in% c("CDK4wt_a", "CDK4wt_b"))) +
geom_point(colour = alpha("darkgreen", 0.7), data = . %>% filter(Id %in% c("CDK4mut_a", "CDK4mut_b"))) +
theme_light(base_size = 11) +
ylim(-5, 1) +
xlim(7, 14) +
xlab(expression(paste("Mean normalized counts ", '(', log[2], ")"))) +
ylab("")
d1000.plot
# Plot as single panel
d10.plot + d100.plot + d300.plot + d1000.plot
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggpubr_0.4.0 forcats_0.5.1
## [3] stringr_1.4.0 dplyr_1.0.5
## [5] purrr_0.3.4 readr_1.4.0
## [7] tidyr_1.1.3 tibble_3.1.0
## [9] ggplot2_3.3.3 tidyverse_1.3.0
## [11] readxl_1.3.1 DESeq2_1.30.1
## [13] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [15] MatrixGenerics_1.2.1 matrixStats_0.58.0
## [17] GenomicRanges_1.42.0 GenomeInfoDb_1.26.4
## [19] IRanges_2.24.1 S4Vectors_0.28.1
## [21] BiocGenerics_0.36.0 ggsci_2.9
## [23] patchwork_1.1.1 knitr_1.31
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 fs_1.5.0 lubridate_1.7.10
## [4] bit64_4.0.5 RColorBrewer_1.1-2 httr_1.4.2
## [7] tools_4.0.3 backports_1.2.1 bslib_0.2.4
## [10] utf8_1.2.1 R6_2.5.0 DBI_1.1.1
## [13] colorspace_2.0-0 withr_2.4.1 tidyselect_1.1.0
## [16] curl_4.3 bit_4.0.4 compiler_4.0.3
## [19] cli_2.3.1 rvest_1.0.0 xml2_1.3.2
## [22] DelayedArray_0.16.3 labeling_0.4.2 sass_0.3.1
## [25] scales_1.1.1 genefilter_1.72.1 digest_0.6.27
## [28] foreign_0.8-80 rmarkdown_2.7 rio_0.5.26
## [31] XVector_0.30.0 pkgconfig_2.0.3 htmltools_0.5.1.1
## [34] highr_0.8 dbplyr_2.1.0 fastmap_1.1.0
## [37] rlang_0.4.10 rstudioapi_0.13 RSQLite_2.2.5
## [40] farver_2.1.0 jquerylib_0.1.3 generics_0.1.0
## [43] jsonlite_1.7.2 BiocParallel_1.24.1 zip_2.1.1
## [46] car_3.0-10 RCurl_1.98-1.3 magrittr_2.0.1
## [49] GenomeInfoDbData_1.2.4 Matrix_1.2-18 Rcpp_1.0.6
## [52] munsell_0.5.0 fansi_0.4.2 abind_1.4-5
## [55] lifecycle_1.0.0 stringi_1.5.3 yaml_2.2.1
## [58] carData_3.0-4 zlibbioc_1.36.0 grid_4.0.3
## [61] blob_1.2.1 crayon_1.4.1 lattice_0.20-41
## [64] haven_2.3.1 splines_4.0.3 annotate_1.68.0
## [67] hms_1.0.0 locfit_1.5-9.4 pillar_1.5.1
## [70] ggsignif_0.6.1 geneplotter_1.68.0 reprex_1.0.0
## [73] XML_3.99-0.6 glue_1.4.2 evaluate_0.14
## [76] data.table_1.14.0 modelr_0.1.8 vctrs_0.3.7
## [79] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [82] cachem_1.0.4 openxlsx_4.2.3 xfun_0.22
## [85] xtable_1.8-4 broom_0.7.5 rstatix_0.7.0
## [88] survival_3.2-7 AnnotationDbi_1.52.0 memoise_2.0.0
## [91] ellipsis_0.3.1